A spatially resolved timeline of the human maternal–fetal interface

Beginning in the first trimester, fetally derived extravillous trophoblasts (EVTs) invade the uterus and remodel its spiral arteries, transforming them into large, dilated blood vessels. Several mechanisms have been proposed to explain how EVTs coordinate with the maternal decidua to promote a tissue microenvironment conducive to spiral artery remodelling (SAR)1–3. However, it remains a matter of debate regarding which immune and stromal cells participate in these interactions and how this evolves with respect to gestational age. Here we used a multiomics approach, combining the strengths of spatial proteomics and transcriptomics, to construct a spatiotemporal atlas of the human maternal–fetal interface in the first half of pregnancy. We used multiplexed ion beam imaging by time-of-flight and a 37-plex antibody panel to analyse around 500,000 cells and 588 arteries within intact decidua from 66 individuals between 6 and 20 weeks of gestation, integrating this dataset with co-registered transcriptomics profiles. Gestational age substantially influenced the frequency of maternal immune and stromal cells, with tolerogenic subsets expressing CD206, CD163, TIM-3, galectin-9 and IDO-1 becoming increasingly enriched and colocalized at later time points. By contrast, SAR progression preferentially correlated with EVT invasion and was transcriptionally defined by 78 gene ontology pathways exhibiting distinct monotonic and biphasic trends. Last, we developed an integrated model of SAR whereby invasion is accompanied by the upregulation of pro-angiogenic, immunoregulatory EVT programmes that promote interactions with the vascular endothelium while avoiding the activation of maternal immune cells.


Supplementary Tables Guide
-Patients and block histology. First tab-Patient meta-data such as age, ethnicity, body mass index, parity and relevant medical conditions such as HIV. Second tab-Information on the histological characteristics of the blocks retrieved, including the presence of cell column anchoring villi. Per patient's block, indication whether cell column anchoring villi were present (1) or absent (0), and the number of regions containing spiral arteries annotated as appropriate for TMA construction by the pathologist (Methods). In blocks containing ≥ 2 distinct, separate pieces of tissue, cell column villi were considered present if they were present on any piece containing pathologist annotations.
Supplementary Table 2 -Artery properties and staging. First tab-Arteries meta-data, including their measured digitized morphological features (see Methods), manual stage and remodeling score δ. Second tab-LDA coefficients for artery morphological features. Per feature LDA coefficients (ld1,ld2) and their Z scored absolute values.

Supplementary Table 3 -Cell-artery and cell-cell spatial enrichment.
First tab-Linear regression results of cell types enrichment around arteries as function of remodeling score δ. Rows: cell types. Columns: regression R 2 , p-value, maximal absolute value of enrichment and trend size (see Methods). Second tab -temporal trends in cell-cell enrichment. Rows: combinations of 2 cell types in the format: cell type a-cell type b. The information in the row is for the enrichment of cell type a around cell type b. The columns show values for the linear regression on per image basis of enrichment scores as function of GA and remodeling score δ: R 2 , the maximal obtained regression R 2 , p-value, maximal enrichment (absolute value) and trend size (see Methods). Third tab-constant in time cell-cell enrichment. Per marker antibody information, including: clone, vendor, vendor ID, channel and elemental reporter, and final staining titers used. The parameters used for marker-specific low-level processing of MIBI data (background removal, denoising, and aggregate removal steps as previously described) are also shown. Second tab-Nanostring morphology markers information. Per marker antibody information including clone, channel and final staining titers used. Table 8 (Fig. 2b). The columns show values for the linear regression on per image proportions as function of GA and remodeling score δ: the log transformed ratio of R 2 , the maximal obtained regression R 2 (maximal between GA and δ), the minimal obtained regression p-value and trend size (see Methods). Second tab-Values plotted in Fig. 5a. Each row represents a combination of a cell type and a functional marker. The columns show values for the linear regression on per image marker positivity rates for the marker-cell type as function of GA and remodeling score δ: the log transformed ratio of R2, the maximal obtained regression R2 (maximal between GA and δ), the minimal obtained regression p-value and trend size (see Methods).

Supplementary
Supplementary Table 11 -NanoString DSP counts and meta-data. First tab-ROI meta-data and raw counts. Meta-data for ROIs. Second tab-mapping Nanostring ROIs to their matching MIBI images. Third tab-raw counts data per ROI. Fourth tab-Differentially expressed genes between interstitial and intravascular EVTs. p-values calculation and multiple hypothesis correction performed by the R Limma package, See methods.

Antibody Validation
Antibody Panel Single Channel Immune Controls  (100 clusters

Analysis of NanoString Whole Transcriptome Data for Decidua ROIs
For this analysis on decidua samples, only genes with background subtracted, normalized counts ≥10 in at least two decidua samples were considered. This resulted in 13104 expressed genes (see Methods "Normalization and scaling of GeoMx counts data").
To estimate the contribution of cell type frequencies in decidua samples to the observed changes in gene expression, cell-type abundancies were estimated in decidua ROIs using the online version of the deconvolution algorithm CIBERSORTx. A custom signature matrix was created using the droplet-based 10x single-cell transcriptome profiles from 30 maternal-fetal tissue samples from Vento-Tormo et al (dataset: https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-6701#). Only cell annotations originating from decidua and placenta were used in creating the signature matrix. This signature matrix was used to impute cell fractions on the Nanostring decidua ROIs (Supplementary Table 13). We then compared the cell type frequencies imputed by CIBERSORTx to the cell type frequencies we measured by MIBI in the decidua masks of the matched MIBI images. We limited this analysis to cell types that are abundant enough in decidua based on our MIBI measurements (frequency ≥1%). In order to match cell type annotations between the two different data sets, we grouped the cell types from the Vento-Tormo et al. (Nature, 2018) dataset and matched them with MIBI as specified in Supplementary Table 13. We then leveraged the deconvolved cell frequencies to predict which genes are trending with gestational age in decidual samples and which trends are better described by tissue composition or SAR. To do that, for each gene, we performed the following analysis: 1. Perform separate polynomial regressions (same as analysis described for artery samples in Methods) of gene expression as a function of gestational age, SAR and each of the six cell type frequencies described in Supplementary Table 13 2. Considering only regression models with p-values ≤ 0.05 and fold change ≥ 2, rank models by their R 2 3. Classify the gene as trending with gestational age, SAR or with one of the cell frequencies based on the model with maximal R 2 4. If there were no regression models with p-values ≤ 0.05 and fold change ≥ 2 for the gene, classify it as non-trending For Coordinated gene expression by pathways in decidua samples we considered all genes for which GA trend was significantly better than SAR trend regardless of tissue composition trends ( Supplementary Information Figure 3). Analysis was performed the same way as described for arteries in Methods section "Coordinated gene expression by pathways in artery".